home *** CD-ROM | disk | FTP | other *** search
/ NetNews Offline 2 / NetNews Offline Volume 2.iso / news / comp / std / c / 117 < prev    next >
Encoding:
Text File  |  1996-08-06  |  3.1 KB  |  122 lines
  1. Path: news.rmii.com!usenet
  2. From: jcoffin@rmii.com (Jerry Coffin)
  3. Newsgroups: comp.std.c
  4. Subject: Re: Help, best way to compare doubles
  5. Date: Tue, 16 Jan 1996 22:44:58 GMT
  6. Organization: TAEUS
  7. Message-ID: <4dh5v8$5sf@natasha.rmii.com>
  8. References: <4d1k09$aqq@mercury.IntNet.net>
  9. NNTP-Posting-Host: slip8163.rmii.com
  10. X-Newsreader: Forte Free Agent 1.0.82
  11.  
  12. jtomich@IntNet.net (Jeff Tomich) wrote:
  13.  
  14. >Having a hard time to figure out a function on how to compare doubles. 
  15. >Any ideas?
  16.  
  17. After a couple of hints in email, and looking at some of the other
  18. posted solutions (and posted problems) I've come up with the following,
  19. which should avoid most of the problems cited so far.  It should work
  20. correctly over a relatively wide range, while hopefully retaining
  21. reasonable speed:
  22.  
  23. /* fcomp.h */
  24. #ifndef FCOMP_H_INCLUDED
  25. #define FCOMP_H_INCLUDED
  26.  
  27. #include <math.h>
  28. #include "dist.h"
  29.  
  30. int fcomp(double a, double b);
  31.  
  32. #endif
  33.  
  34. /* fcomp.c */
  35. #include "fcomp.h"
  36.  
  37. /* This asks for equality to roughly 10 digits.
  38.  */
  39. static const double delta = 1e-20;
  40.  
  41.  
  42. /* Iterations   Accuracy
  43.  *  2          6.5 digits
  44.  *  3           20 digits
  45.  *  4           62 digits
  46.  * assuming a numeric type able to maintain that degree of accuracy in
  47.  * the individual operations.
  48.  */
  49. #define ITER 3
  50.  
  51. static double dist(double P, double Q) {
  52. /* A reasonably robust method of calculating `sqrt(P*P + Q*Q)'
  53.  *
  54.  * Transliterated from _More Programming Pearls, Confessions of a Coder_
  55.  * by Jon Bentley, pg. 156. (copr. Bell Telephone Labs, 1988, published
  56.  * Addison-Wesley.  ISBN 0-201-11889-0 )
  57.  */
  58.  
  59.     double R;
  60.     int i;
  61.  
  62.     P = fabs(P);
  63.     Q = fabs(Q);
  64.  
  65.     if (P<Q) {
  66.         R = P;
  67.         P = Q;
  68.         Q = R;
  69.     }
  70.  
  71. /* The book has this as:
  72.  *  if P = 0.0 return Q; # in AWK
  73.  * However, this makes no sense to me - we've just insured that P>=Q, so
  74.  * P==0 only if Q==0;  OTOH, if Q==0, then distance == P...
  75.  */
  76.     if ( Q == 0.0 )
  77.         return P;
  78.  
  79.     for (i=0;i<ITER;i++) {
  80.         R = Q / P;
  81.         R + R * R;
  82.         R = R / (4.0 + R);
  83.         P = P + 2.0 * R * P;
  84.         Q = Q * R;
  85.     }
  86.     return P;
  87. }
  88.  
  89. int fcomp(double a, double b) {
  90. /* Hopefully a reasonably robust comparison of a and b.
  91.  * Returns 1 if they are equal to the precision indicated by delta.
  92.  * Returns 0 if the two differ by more than delta.
  93.  * As far as I can tell, the major possibility of overflow is if a and b
  94.  * are both within a factor of 2 of DBL_MAX in magnitude, and differ in
  95.  * sign.  
  96.  * Delta is small enought that b+delta should never overflow: if  b is 
  97.  * is close to DBL_MAX, adding delta will generally have no effect.
  98.  */
  99.  
  100.     double difference = fabs(a-b);
  101.  
  102.     double distance = dist(a,b+delta);
  103.  
  104.     if ( difference / distance < delta )
  105.         return 1;
  106.     return 0;
  107. }
  108.  
  109. I'm _NOT_ particular skilled at numeric programming, so anybody who is
  110. is more than welcome to tear it up if they see fit - this is basically
  111. just one of the other posts recast to use Moler & Morrison's distance
  112. finding algorithm...
  113.     Later,
  114.     Jerry.
  115.  
  116. /* I can barely express my own opinions; I certainly can't
  117.  * express anybody else's.
  118.  *
  119.  * The universe is a figment of its own imagination.
  120.  */
  121.  
  122.